Exotic spin orders driven by orbital fluctuations in the Kugel-Khomskii model 
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We study zero temperature phase diagram of the three-dimensional Kugel-Khomskii model on a 
cubic lattice using the cluster mean field theory and different perturbative expansions in the orbital 
sector. The phase diagram is rich, goes beyond the single-site mean field theory due to spin-orbital 
entanglement. In addition to the antiferromagnetic (AF) and ferromagnetic (FM) phases, one finds 
also a plaquette valence-bond phase with singlets ordered either on horizontal or vertical bonds. 
More importantly, for increasing Hund's exchange we identify three phases with exotic magnetic 
order stabilized by orbital fluctuations in between the AF and FM order: (i) an AF phase with two 
mutually orthogonal antiferromagnets on two sublattices in each ab plane and AF order along the c 
axis (ortho-G-type phase), (ii) a canted- ^4-type AF phase with a non-trivial canting angle between 
nearest neighbor FM layers along the c axis, and (iii) a striped- AF phase with anisotropic AF order 
in the ab planes. We elucidate the mechanism responsible for each of the above phases by deriving 
effective spin models which involve second and third neighbor Heisenberg interactions as well as 
four-site spin interactions going beyond Heisenberg physics, and explain how the entangled nearest 
neighbor spin-orbital superexchange generates spin interactions between more distant spins. 
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I. INTRODUCTION 

Recent interest and progress in the theory of spin- 
orbital superexchange models was triggered by the obser- 
vation that orbital degeneracy drastically increases quan- 
tum fluctuations (QF) which may suppress long-range or- 
der in the regime of strong competition between different 
types of ordered states near the quantum critical pointi 
The simplest and archetypical three-dimensional (3D) 
model for the spin-orbital physics is the Kugel-Khomskii 
(KK) model introduced for KCUF3 by Kugel and Khom- 
skii long ago<2 KCUF3 is a strongly correlated system 
with a single hole within degenerate e g orbitals at each 
Cu 2+ ion. Kugel and Khomskii showed that interactions 
between strongly correlated electrons could then stabi- 
lize orbital order by a purely electronic mechanism, the 
superexchange. A similar situation occurs in a number of 
compounds with active orbital degrees of freedom, where 
strong on-site Coulomb interactions localize electrons (or 
holes) and give rise to spin-orbital superexchange^— 

The orbital superexchange may stabilize the orbital 
order by itself, but in e g systems it is usually helped 
by the Jahn- Teller distortions of the lattice which gen- 
erate effective intersite orbital interactions For in- 
stance, in LaMnOs the terms that originate from the su- 
perexchange and the Jahn- Teller distortions are of equal 
importance and both are necessary to explain the ob- 
served high temperature Too — 780K of the structural 
transition^ Also in KCUF3 the lattice distortions play an 
important roleii~— and are responsible for its strongly 
anisotropic magnetic and optical properties. The latest 
theoretical and experimental results for this compound 
show that other types of interactions, like Goodenough 
processes in the superexchange^ direct orbital exchange 



driven by a combination of electron-electron interactions 
and ligand distortions^ or dynamical Dzyaloshinsky- 
Moriya interaction^ are necessary to explain structural 
phase transition in KCUF3 at Too — 800.fr and peculiar- 
ities in spin dynamics^ The compound itself is believed 
to be the best realization of the one-dimensional (ID) an- 
tiferromagnetic (AF) Heisenberg model above the Neel 
temperature Tv = 39K^ and spinon excitations were 
indeed observed in neutron scattering^ 

While the coexisting A-type AF (A-AF) order and the 
orbital order are well established in KCUF3 below the 
Neel temperature Tjy — 39 K^ and this phase is re- 
produced by the spin-orbital d 9 superexchange model in 
the mean field (MF) approximation^ the phase diagram 
of this model is still unknown beyond the MF approach 
because of strongly coupled spin and orbital degrees of 
freedomii 2 ^ which poses an outstanding question in the 
theory: Which types of coexisting spin and orbital order 
(or disorder) are possible when its microscopic parame- 
ters (splitting of the e g orbitals E z and Hund's exchange 
Jh) arc varied? So far, it was suggested that the long- 
range AF order is destroyed by spin-orbital QF,— >2£ but 
another possibility that the ordered state might be sta- 
bilized by the order-out-of-disorder mechanism was also 
pointed out in the regime of ferro-orbital (FO) order 
of 2>z 2 — r 2 orbitals^ An alternative to magnetic order 
are spin-disordered phases with pronounced valence-bond 
spin-orbital correlations, as suggested by simple varia- 
tional wave functions X 

The purpose of this paper is to investigate the phase 
diagram of the 3D KK model, in the presence of spin- 
orbital QF. This subject is of interest in the broad con- 
text of frustration in magnetic systems ; 22 ! 23 which ap- 
pears to be intrinsic for the orbital superexchange 
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To establish reliable results concerning short-range or- 
der in the crossover regime between phases with long- 
range AF or ferromagnetic (FM) order, we developed a 
cluster MF approach which goes beyond the single-site 
MF in the spin-orbital system 2 ^ and is based on an ex- 
act diagonalization of a cluster coupled to its neighbors 
by the MF terms. The cluster is chosen to be suffi- 
cient for investigating both AF phases with four sub- 
lattices and valence-bond states, with spin singlets ei- 
ther along the c axis or within the ah planes. This the- 
oretical approach is motivated by possible spin-orbital 
entanglement^^ which is particularly pronounced in the 
ID SU(4) [or SU(2)®SU(2)] spin-orbital model 2 ! and oc- 
curs also on the frustrated triangular lattice,— and in the 
models for perovskites when spin correlations are AF on 
the bonds! — then the Goodcnough-Kanamori rules^ 
are violated in some cases. In the perovskite vanadates 
such entangled states play an important role at finite 
temperature: in their optical properties^ in the phase 
diagram^! and in the dimerization of FM interactions 
along the c axis in the C-AF phase of YV03^ under- 
stood within the ID spin-orbital models Previous stud- 
ies employing the cluster MF approach have shown that 
phases with entangled spin-orbital degrees of freedom oc- 
cur in the KK bilayer^ while a noncollinear spin order 
emerges from entangled spin-orbital fluctuations in the 
two-dimensional (2D) KK monolayer^ Below we inves- 
tigate whether spin-orbital entangled states could also 
play a role in the present 3D KK model. 

The paper is organized as follows. In Sec. Ill Al we 
briefly introduce the spin-orbital superexchange in the 
KK model. As the first approximation, in Sec. Ill Bl we 
present the phase diagram of the model obtained in the 
single-site MF approximation — this approach ignores 
any possible spin-orbital entanglement. To include the 
effects of spin-orbital QF we introduce the cluster MF 
method for the 3D system in Sec. IIII Al with two differ- 
ent topologies of the cluster, and present the phase dia- 
gram modified by spin-orbital entanglement in Sec. IIII Al 
It contains three phases with exotic magnetic order: the 
ortho-G-AF phase similar to that found recently for the 
monolayer^ the canted- A- AF phase, and the striped- AF 
phase. Next we present the behavior of the order param- 
eters, spin angles, total magnetization, correlations and 
spin-orbital covariances for two paths in the phase dia- 
gram which include certain exotic phases: (i) from the 
A-AF through canted- A-AF to FM phase in Sec. iHTCl 
and (ii) from the striped- AF to G-AF phase in Sec. IIII Dl 
As we show in Sec. IIII El rather weak orbital fluctuations 
are found in several phases and the orbital moments are 
reduced by them. Then we explain and highlight the ori- 
gin of the exotic magnetic phases by deriving effective 
spin Hamiltonians within perturbative expansion in the 
orbital sector. Spin models for the ortho-G-AF phase, 
canted- A-AF phase, and striped- AF phase are derived in 
Sees. ITVAl IPTBl and HVDl respectively In Sec. HVCl 
using a similar perturbative expansion, we explain the 
absence of the G-AF phase in the phase diagram of the 



KK model obtained within the cluster MF approxima- 
tion. Summary and conclusions are given in Sec. |Vl 
while certain additional details of the performed pertur- 
bative analysis are presented in Appendices A-C. 

II. THE KUGEL— KHOMSKII MODEL 
A. Frustrated spin-orbital superexchange 

For realistic parameters the late transition metal ox- 
ides or fluorides are strongly correlated and electrons lo- 
calize in the 3d orbitals t 36 i 37 leading in cuprates to Cu 2+ 
ions with spin S — 1/2 in d 9 configuration and e g orbitals 
occupied by one hole: 

\x) = (x 2 - y 2 )/V2, \z) = [3z 2 - r 2 )/V6 . (2.1) 

The examples of such systems are: KCUF3 with 3D cubic 
lattice, K3CU2F7 representing bilayer compounds, and 
K2CUF4 and La2Cu04 with 2D square lattice. We also 
use below the short-hand notation for the orbital basis, 
{x,z}. Here t^g orbitals are split in the octahedral field 
and do not couple to e 9 's by hopping through fluorine, so 
they can be neglected. In what follows we investigate the 
3D spin-orbital superexchange model obtained by consid- 
ering charge excitations between transition metal ions^ 
d\ d 9 ^ df dj° in the regime of large U, and neglect the 
coupling to lattice distortions arising due to the Jahn- 
Teller lattice distortions. 

One finds the Heisenberg Hamiltonian for S = | spins 
coupled to the orbital problem, 

n = -\ J E E flS-^E'f. (2 - 2) 

7— a,b,c (ij) 1 17 i 

where bond-interaction terms H?j are defined as follows, 




Here 7 = a, b, c labels the direction of a bond (ij) in the 
3D system. The energy scale is given by the superex- 
change constant, 



where the hopping t is the effective intersite (dda) hop- 
ping element for an e g hole between z orbitals along the 
c axis^ and the other e g hopping elements obey the 
cubic symmetry^ The orbital operators at site i are 
Ti = {t°-,t^, rf }. These operators represent e g orbital de- 
grees of freedom, have cubic symmetry in the 3D lattice 
and are expressed in terms of Pauli matrices {erf, erf, of} 
in the following way: 

rf h) = l(-af±V3af), Tf = ±of. (2.5) 
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The matrices {cr 7 } act in the orbital space (and have 
nothing to do with the physical spin S^ present in this 
problem). Note that {t?} operators are not indepen- 
dent from one another because they satisfy the local con- 
straint, J2j T i = 0- The operators IT^ and H'^ stand for 
projections of spin states on the bond (ij) on a singlet 
(Ilf ■) and triplet (II' ) configuration, respectively, 



n 



- s, • s. 



n 



(2.6) 

for spins S = h at both sites i and j of the bond (ij). 
Spin interactions obey the SU(2) symmetry. More details 
on the derivation of Eq. (|2.3[) may be found in Ref. [l^. 

The model Eq. (|2.2[) depends thus on two 
parameters*^ (i) Hund's exchange coupling 77, and (ii) 
the crystal- field splitting of e g orbitals E z /J. The first 
of them is given by the ratio of Hund's exchange Jh and 
intraorbital Coulomb element U defined in a standard 
way as in the degenerate Hubbard model^S 



Jh 
U ' 



(2.7) 




FIG. 1. (Color online) Left panel: schematic view of four 
representative orbital configurations on a representative cube 
of the 3D lattice: (a) AO order with (t"' 6 ') = 1/2 changing 
from site to site and (rf) = —1/4, obtained for E z < 0, (b) 
AO order with (t"' 6 ') = —1/2 changing from site to site and 
(rf) — 1/4, obtained for E z > 0, (c) FO order with occupied z 
orbitals and (rf) = —1/2 (cigar-shaped orbitals), and (d) FO 
order with occupied x orbitals and (rf) = 1/2 (clover-shaped 
orbitals). Right panel: schematic view of four representative 
spin configurations (arrows stand for up or down spins) shown 
on a cubic cluster: (i) A-AF configuration, (ii) C-AF config- 
uration, (iii) FM configuration, and (iv) G-AF configuration. 



and determines the values of the coefficients in Eq. (|2.3 
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1-77' 
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V 



(2.8) 



The typical energies for the Coulomb U and Hund's ex- 
change Jh can be deduced from the atomic spectra or 
from density functional theory with constrained electron 
densities. Earlier studies performed within the local den- 
sity approximation (LDA) gave rather large values of the 
interaction parameters for Cu 2+ ions>2I U = 8.96 eV and 
Jh = 1.19 cV. More recent studies used the Coulomb in- 
teractions treated within the LDA+U scheme and gave 
somewhat reduced values^ - U = 7.5 eV and Jh = 0.9 
eV. However, both parameter sets give a rather similar 
value of Hund's exchange parameter 77, being within the 
expected range 0.10 < rj < 0.15 for strongly correlated 
late transition metal oxides.— Note that the physically 
acceptable range is much broader, i.e., < 77 < 1/3. The 
upper limit follows from the condition (U — 3Jh) > for 
the high-spin excitation energy. 

The last term of the H e Hamiltonian (|2.2p lifts the 
degeneracy of the two e g orbitals, 



H 7 



(2.9) 



where \ni xai rii za \ are hole number operators in 
orbitals (|2.1I) at site i. It favors hole occupancy of a; 
(z) orbitals when E z > (E z < 0) and can be associated 
with a uniaxial pressure along the c axis or a crystal-field 
splitting induced by a static Jahn- Teller effect. 

In Figs. QJa)-QJd) we present typical orbital configu- 
rations with FO order and alternating orbital (AO) or- 
der considered in the e g orbital models! 10 ' 38 In the next 
sections we analyze their possible coexistence with spin 



order in the KK model Eq. (|2.2p . As we can see, the 
maximal (minimal) value of the orbital operators is 
related with orbital taking shape of a clover (cigar) with 
symmetry axis pointing along the direction 7. 



B. Phase diagram in single-site mean field 

After averaging over spins, the Hamiltonian of Eq. 
(|2.2j) . originally expressed in terms of bond operators, 
can be rewritten as an effective orbital Hamiltonian, 



H 



MF : 



r7 +7 (x 7 -0 + T-7r-^X 7 +0 



(2.10) 



where i runs over sites of the cubic lattice, and i + j is the 
nearest neighbor (NN) of site i along the axis 7 = a, b, c. 
The coefficients, 



X 



nil] + r 2 n] 



C = (r 2 +n)W s , (2.11) 



are parameters obtained by averaging of the spin projec- 
tors in Eq. (|2.6[) under assumption that the spin order 
depends only on the direction 7 and all the bonds (i, i+7) 
along the axis 7 are equivalent: 



H2 = i -(S l .S J+7 ), 



n 



7 _ 



(Si • S ? ; 



+7/ 



(2.12) 



In the single-site MF the spin QF are absent at zero tem- 
perature and the projectors can be replaced by their av- 
erage values. This is sufficient to investigate the phases 
with either AF or FM long-range order. 
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TABLE I. Mean values of triplet 117 an d singlet pro- 
jection operators Eqs. (|2.12[) for a bond (ij) along the axis 
7 = a, b, c in different phases with long-range magnetic order 
which occur in the MF phase diagram, see Figs. [Hi)-[]Jiv). 
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G-AF 


1/2 


1/2 


1/2 


1/2 


C-AF 


1/2 


1 


1/2 





A-AF 


1 


1/2 





1/2 


FM 


1 


1 









The values of the projection operators (|2.12j) depend 
on the assumed spin order. Here we consider four differ- 
ent spin configurations shown in Fig. Q] (i) A-AF phase 
— with FM order in the ab planes and AF correlations 
along the c axis [Fig. H^i)], (ii) C-AF phase — with AF 
order in the ab planes and FM correlations along the c 
axis [Fig. [I]n)], (iii) FM phase [Fig. [TJiii)], and (iv) G- 
AF phase Neel state [Fig. [T[iv)]. In the single-site MF 
approximation we use the classical average values of the 
spin projection operators in the above phases listed in 
Tabic HI Apart from fully AF and FM phase we include 
also A-AF (C-AF) configurations with spin correlations 
being AF along the c axis (in the ab planes) , and FM oth- 
erwise. Solutions of the self-consistency equations and 
ground state energies in different phases can be obtained 
analytically, as shown in Ref. HH- 

The phase diagram presented in Fig. [5] was obtained 
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z 

FIG. 2. (Color online) Phase diagram of the 3D KK model 
obtained in the single-site MF approximation. Shaded dark 
gray (green) area indicates phases with AO order while the 
remaining magnetic phases are accompanied by FO order with 
fully polarized orbitals, either x (for E z > 0) or z (for E z < 0). 
In this approach the G-AFx and C-AFa; phases with FO order 
are degenerate. 



by purely energetic considerations — it shows the bor- 
der lines between phases with the lowest energies in the 
(E z , 77) plane. Remarkably, this phase diagram is almost 
the same as the one obtained for the bilayer system^ 

- the same phases were found both in the 3D and 
in the bilayer KK model, and the major difference is 
the location of the multicritical point, being now at 
(E z ,rj) = (0,0). This reflects the cubic symmetry in 
the model (|2.2|) at E z — 0, while this symmetry is bro- 
ken by the crystal-field term in both planar models - 
indeed, for the KK bilayer the multicritical point is lo- 
cated at (E Zl rj) = (— 0.25 J, 0)j24 and it is moved further 
to (E z ,Tj) = (— 0.5 J, 0) for the monolayer KK models 
At ?; = one finds only two AF phases: (i) G-AFz for 
E z < and (ii) G-AFx for E z > 0, both with polar- 
ized orbital configuration (FO order) which involves ei- 
ther cigar-shaped z orbitals in the G-AFz phase, see Fig. 
QJc), or clover-shaped x orbitals in the G-AFx, see Fig. 
[Tfd). Because of the planar orbital configuration in the 
G-AFx phase one finds no interplane exchange coupling 

- thus the spin order along the c axis is undetermined 
and this phase is degenerate with the C-AF one. This 
degeneracy is lifted in the cluster MF approach, see be- 
low. 

For higher r\ the number of phases increases abruptly 
by those with AO configurations (green areas), as shown 
in Figs. [TJa) and[IJb), coexisting with different possible 
spin orders: the A-AF, C-AF, G-AF and FM phases, re- 
spectively. Altogether, the phase diagram obtained here 
reproduces qualitatively the one obtained before for the 
simplified model using the lowest order expansion in r\X 
The C-AF phase occurs in a narrow range of parameters 
in between the G-AF and either A-AF or FM phase. 

In contrast to the FO phases, the AO order in the 
shaded phases is never trivial in the sense that the or- 
bitals are never fully polarized in any direction, which 
is a feature of the self-consistent MF solution (see Ref. 
[HI). The only new phase with FO configuration is the 
A-AFz one appearing above 77 = 0.155 for E z < 0. On 
the other hand, the FM spin order coexists solely with 
alternating orbitals. Finally, the new exotic magnetic 
phases reported in Sec. IIIII (ortho-G-A, canted- A-AF, 
and striped- AF phase) cannot appear here as their sta- 
bilizing mechanism is absent in the single-site MF, so 
they are not included in Table [IJ 



III. CLUSTER MEAN FIELD APPROACH 

A. Clusters and order parameters 

Following the ideas from Ref. Hi| and [3^, to obtain 
more insight into the phase diagram of the 3D KK model 
and to include the QF and spin-orbital entanglement on 
the bonds, we divide the lattice into clusters, either pla- 
quettes shown in Fig. EJa), or chain clusters shown in 
Fig. GUb). The most natural choice of the cluster would 
be a cube with eight sites as it was done in the bilayer 
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easels but, since we want to keep two components of the 
spin order parameter, we adopted here a simpler and less 
time-consuming approach using four-site clusters. The 
geometry of the selected clusters depends on the direc- 
tion in which we expect a large energy gain due to QF. 
For example, if we want to study the transition between 
the G-AF and A-AF phases then the reasonable cluster 
topology is a square in the ab plane, shown in Fig. EJa), 
because the order along the c axis does not change across 
the transition. On the other hand, if we are interested in 
a transition between the A-AF and FM phase where the 
order in the ab planes remains constant, then a better 
choice is a chain along the c axis — see Fig. [3th). 

The interactions along bonds corresponding to the 
solid lines in Fig. [4] are treated by exact diagonalization 
as they stand in the 3D KK model (|2.2j) . while the bonds 
represented by dashed lines are decoupled in the MF ap- 
proximation using the approximate identity for any bond 
operator: 

OiOjKOiW-^OiHOi) 

+ O j (O i )-±(O i )(O j ). (3.1) 

To simulate infinite 3D lattice we need MF bonds in all 
three directions. Now we assume that site i belongs to a 
chosen cluster and j belongs to a neighboring one, and 
the bond is splitted into two halves — the first one is 
added to the Hamiltonian of the cluster i and the other 
one to the cluster j. In this way the original KK Hamil- 
tonian transforms into the sum of commuting cluster 
Hamiltonians interacting via MF terms. 

The MFs follow from Eq. (f3~Tj) applied to all the two- 
site operator products encountered in the Hamiltonian 
(12.21) and are defined as follows: 

s? = {St) , tl = (t7) , << 7 = (S?t7) . (3.2) 

Here a = x, z, 7 = a, b (7 = c) and i = 1, 2, 3, 4 for the 
cluster sites of a plaquette (chain) cluster, see Fig. [3] 
Note that the SU(2) symmetry of the spin sector does 
not need to be broken in the z spin direction (a = z), 
but we also allow a — x to capture more exotic types of 
magnetic order suggested by the results reported recently 
for the 2D system^ However, we do not need to consider 
a = y because the KK Hamiltonian is real. 

To obtain the unbiased and the most general phase 
diagram, we do not assume anything about the order in- 
side the cluster because the considered plaquette is small 
enough to keep the order parameters at all its sites as 
independent variables along the MF iteration process. 
Nevertheless, we still need to relate different clusters to 
one another to make the problem solvable. In the case 
of a square cluster we assume that the neighboring clus- 
ters in the c direction can have either the same spin or 
inverted spin configuration (it gives either FM or AF 
bonds along c axis). Furthermore, we assume that the 
neighbors in the ab plane can have either the same or- 
bital configuration that gives AO and FO orders in the 
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FIG. 3. Schematic view of the clusters (solid lines) used in 
the cluster MF approach of Sec. IHII to the 3D KK model: (a) 
a plaquette in the ab plane, and (b) a chain along the c axis. 
Vertices i = 1,2,3,4 and directions 7 = a,b,c are marked; 
dashed lines stand for the outgoing bonds where the spin- 
orbital interactions are replaced by the MF terms at neigh- 
boring sites, see Eq. (|3.ip . 



ab planes, or the orbital configuration is rotated by 7r/2 
in the ab plane — it gives the plaquette valence-bond 
(PVB) phase. Similarly, we assume that the chain clus- 
ters are copied without any change along the c axis and 
the neighboring chains in the ab planes have: (i) orbital 
configuration rotated by 7r/2, and (ii) spin configuration 
either inverted (it gives planar AF order) or unchanged 
(it gives planar FM order). All these assumptions are 
necessary to solve the self-consistent cluster MF problem 
and are motivated by the phase diagrams of the bilayer 
and the monolayer systems (see Refs. and[H). 

The self-consistency equations still cannot be solved 
exactly because the effective cluster Hilbert space is of 
the size d = 2 8 which is too large for analytical methods. 
The way out is to use Bethe-Peierls- Weiss method, i.e., 
to set certain initial values for the order parameters Eq. 
(|3.2[) , {sf, tj, vf' 7 }, and next to employ numerical di- 
agonalization algorithm to the cluster Hamiltonian. We 
recalculate the order parameters, {sf , t] , vf' 7 }, along the 
iteration process, and this procedure is repeated until the 
convergence conditions for energy and order parameters 
are satisfied. 



B. Phase diagram 

The phase diagram obtained in the cluster MF ap- 
proach is shown in Fig. @] One finds the phases with 
magnetic long-range order, obtained in the single-site 
MF and explained in Sec. Ill Fit in the broad (unshaded) 
part of the phase diagram: the G-AF, A-AF and FM 
phase. The shading in the center marks the spin disor- 
dered PVB phase with pairs of spin singlets alternating in 
the ab planes and accompanied by z-like 3x 2 — r 2 /3y 2 — r 2 
orbitals pointing along the singlet bonds.— Analogous 
valence-bond phases were also found in the bilayer— and 
monolayer— KK model. The darker (orange) shading in- 
dicates exotic magnetic orders which can be found when 
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FIG. 5. (Color online) Classical view of the ortho-G-AF spin 
order realized in the 3D KK model within ab planes, and 
staggered along the c axis. Vertical and horizontal arrows 
correspond to two interpenetrating AF states on the sublat- 
tices of next nearest neighbor sites. Up (down) arrows stand 
for (Sf) = ±1/2, right (left) arrows stand for (Sf) = ±1/2. 



FIG. 4. (Color online) Phase diagram of the 3D KK model 
in the cluster MF approximation. Plaquette valence-bond 
(PVB) phase with alternating spin singlets in the ab planes, 
highlighted in light gray (yellow), occurs between the phases 
with magnetic long-range order, see Fig. [1] Phases with 
exotic magnetic order are shaded in dark gray (orange). 



some AF spin interactions change into FM ones; they are: 
(i) ortho-G-AF phase already encountered in the 2D KK 
model and called there ortho-AF phase,— (ii) canted-A- 
AF phase, and (iii) striped-AF phase. These new phases 
arise from orbital fluctuations in the regimes of strongly 
frustrated spin-orbital superexchange, as explained be- 
low. 

We begin with the ortho-G-AF phase, with the spin or- 
der consisting of two interpenetrating AF sublattices in 
the ab planes, as shown in Fig. [5l This configuration re- 
peats itself in the next ab plane but all spins are inverted, 
meaning AF order along the c axis. Note that this phase 
separates phases with antiferromagnetism (G-AF) and 
ferromagnetism (A-AF) within ab planes, as found be- 
fore in the 2D KK model^ The interactions along the c 
axis are compatible with the in-plane magnetic order and 
even stabilize it as one finds here the ortho-G-AF phase 
at a given r\ for a lower value of E z than that in the 2D 
phase. The ortho-G-AF phase with interplanar antifer- 
romagnetism replaces here the resonating valence-bond 
(RVB) phase found before in the 3D KK model,— where 
it was proposed as an intermediate phase separating the 
G-AF and A-AF phases. We believe that the present re- 
sult is more realistic (at zero temperature) than the RVB 
phase within ID chains along the c axis found before^ — 
this latter phase would be easily modified by any in-plane 
magnetic order because the ID Heisenberg antifcrromag- 
net is critical and thus easily destabilized. In addition, 
orbital fluctuations remove locally AF spin coupling and 
thus block the resonance in the RVB phase. We argue 



below that the ortho-G-AF phase is well justified by the 
effective perturbative spin model derived for the 2D KK 
model in Rcf. [35| and for the present 3D model in Sec. 

As r\ is further increased in the A-AF phase, one finds a 
second magnetic transition, with spin correlations along 
the c axis changing sign. Here the canted-A-AF phase is 
found as an intermediate phase connecting smoothly (in 
contrast to the ortho-G-AF) the A-AF phase with the 
FM one. In the canted-A-AF configuration the spins are 
FM in the ab planes and the order along the direction c 
changes gradually from AF to FM with intcrplane spin 
angle 9, being the canting angle and taking values be- 
tween 6 = and 8 = tt, see Fig. [6jb). On the other 
side of the phase diagram, i.e., for E z > 0, one finds 
the striped-AF phase characterized by symmetry break- 
ing between the a and b directions in the orbital and 
spin sectors, for similar values of rj. The magnetic order 
in striped-AF phase is AF with anisotropy; along one 
direction in the ab plane the order is purely AF and in 
the perpendicular direction the angle between neighbor- 
ing spins is close to (but not exactly) n as shown in Fig. 
[SJa). The orbital configuration is FO with one preferred 
direction, i.e., t a 7^ t b . Striped-AF phase connects with 
left G-AF phase by a smooth phase transition. Further 
on we will present some analytical arguments explaining 
both canted-A-AF and striped-AF phase by perturbative 
expansion, see Sec. II V Bl and Appendix B and Sec. IIVDI 
and Appendix C, respectively. 

Otherwise, the phase diagram of Fig. @] contains the 
G-AF, A-AF and FM configurations placed similarly as 
in the single-site MF phase diagram of Fig. O The de- 
generacy between the left G-AF and G-AF phase is now 
removed, as in the bilayer KK model^i and this time 
we provide a perturbative explanation of this fact in Sec. 
IIV CI Similarly to the bilayer phase diagram, the PVB 
phase connects with the left G-AF phase by the inter- 
mediate PVB-AF configuration, but due to the presence 
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\ 



(b) 



FIG. 6. (Color online) Schematic views of the two exotic spin 
orders realized by the 3D KK model at large Hund's exchange 
r\ > 0.2: (a) striped-AF order in the ab plane, with AF order 
along the c axis and angle (f> a between the NN spins along the 
a axis; (b) spin order realized in the canted- A-AF phase with 
FM order in ab planes and spin canting angle 8 along the c 
axis. The regions of stability of these phases are shown in 
Fig. |4]by dark gray (orange) shading. 



of the right G-AF phase we have also the right PVB- 
AF phase. One can summarize that the whole bottom 
part of the phase diagram up to r\ w 0.085 contains only 
smooth (second order) phase transitions when E z is var- 
ied. Before deriving the effective spin models for the new 
magnetic configurations found in the 3D KK model we 
will look more closely at the phase transitions along two 
cuts in the phase diagram of Fig. @] (i) connecting the 
A-AF and FM phases through the canted- A-AF phase 
(Sec. IOC), and (ii) from the striped-AF to the G-AF 
phase (Sec. HITdI) . 



C. From the A-AF to FM phase 

First, we consider the negative crystal-field splitting 
E z = —0.5 J — for this representative value the order 
changes first from the A-AF into the canted- A-AF phase, 
and next into the FM phase when Hund's exchange rj 
increases. We selected the chain cluster of Fig. EJb) to 
study these phase transitions as the spin order in the ab 
planes does not change. The changes of spin order along 
the c axis are captured by the cosine of the spin canting 
angle 9 along the c axis and the total magnetization |s|, 
defined in the following way: 



cosfl = \ (s\s% + sfsf), 
\s\ = VW + (* Z ) 2 , 



(3.3) 
(3.4) 



and displayed in Figs. [TJa) andJTJb). For the AF configu- 
ration (in the A-AF phase) one finds 9 = it (cos 9 = — 1), 
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FIG. 7. (Color online) Evolution of the magnetic order for 
increasing r\ at E z — — 0.5 J: (a) cosine of the canting angle 



9 (|3.3[) . and (b) total magnetization \s\ (|3.4|) : both in the A- 
AF, canted- A-AF and FM phases from left to right. Quantum 
phase transitions are indicated by vertical dotted lines. 



while for the FM order 9 = (cos 9 = 1). In the canted- 
A-AF phase cos 9 interpolates smoothly between these 
two limiting values. Figure[7Jb) shows that the spin order 
parameter \s\ is gradually reduced and the QF increase 
when rj decreases and the A-AF phase is approached, but 
even in the A-AF phase the spin order is almost classi- 
cal with \s\ ~ 0.5. Indeed, the QF in the A-AF phase 
with all the bonds in ab planes being FM are expected 
to be considerably reduced from the 2D Hciscnbcrg anti- 
ferromagnet, as shown in the spin- wave theory— In the 
canted- A-AF phase the slope of \s\ is the largest and the 
QF almost saturate when the spins have rotated com- 
pletely to the A-AF phase. 
The spin correlations, 



G s c (d) = (S 1 -S 1+d ) 



(3.5) 



along the c axis for the NN, NNN and 3NN (distance 
d = 1, 2, 3) are shown in Fig. [5]^ a), for the same path in 
the parameter space as in Fig. \7\ The NN and 3NN cor- 
relations confirm that the order along the c axis changes 
from the AF to FM one in a continuous way, with spin 
correlations passing through zero. The NNN spin cor- 
relation stays FM and is almost constant in the entire 
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phase. Quantum phase transitions are indicated by vertical 
dotted lines. 
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parameters, and (b) total magnetization \s\ (|3.4[) and cosines 
of spin angles (f> a and 4> b y see Fig- Oa). Quantum phase 
transition is indicated by vertical dotted lines. 



range of rj. This peculiar behavior will be explained by 
an effective perturbative spin Hamiltonian in Sec. IIVBI 
Figure [SJ^b) presents the spin-orbital covariances: the 
on-site ones, 

rl=v^-s*tl, (3.6) 

considered here only for the z spin component sf, and 
the bond covariances, 

R c (d) = ((Sx-Sx+dKr^) - (S 1 -S 1+d ) (r^ 1+d ) , 

(3.7) 

for the NN (d = 1) and for further neighbor (d = 2, 3) 
operators in the chain cluster of Fig. [H^b). As one ex- 
pects, in the FM phase all the covariances vanish and the 
factorization of spin and orbital operators is exact. 

In the canted-A-AF phase the slopes of the covariances 
are the steepest. The R c (l) function is the one of the 
largest magnitude meaning that the spin-orbital entan- 
glement on the NN bonds^ is high. In contrast to that, 
the covariances R c (d) are close to zero for d > 1, but 
i? c (3) as the only one becomes more significant in the 



canted-A-AF phase, indicating that this phase can be 
governed by longer range spin interactions accompanied 
by orbital fluctuations. The on-site covariances behave 
in a monotonous way and reach relatively small absolute 
values meaning that they are not of the prime importance 
in the considered phases. 



D. From the striped-AF to G-AF phase 

Another exotic type of magnetic order found in the 
3D KK model is the striped-AF phase. This phase can 
evolve smoothly towards the ordinary G-AF Neel order 
when E z increases. Here we use the plaquette cluster of 
Fig. [3{a) as the spin order in the ab planes changes. In 
Fig. EI 3 *) we present the evolution of the order parame- 
ters near this transition at r\ — 0.22. The orbital order 
parameters {t a , t b } confirm breaking of the a-b symmetry 
in the striped-AF phase where they take slightly different 
values; this difference vanishes at the phase transition. In 
the magnetic sector we can distinguish four spin sublat- 
tices, see Fig. HJa), two of which are not related by a 
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FIG. 10. (Color online) Evolution of covariances from the 
striped-AF to G-AF phase at r\ — 0.22 under increasing E z /J 
in the ab planes: (a) on-site spin-orbital covariances (|3.6|) for 
the x spin component, and r®', 3 y and (b) on-site spin- 

orbital covariances (|3.6[1 for the x spin component, rf;°i and 



1(3) 

3.6[) . together with bond spin-orbital covariances 7?° 



'1(3) 

(|3.7p . Quantum phase transitions are indicated by vertical 
dotted lines. 



spin inversion. To show the full complexity of the spin 
order we present spin averages on sites i = 1, 3 and both 
spin components a = x, z. The behavior of curves con- 
firms the striped character of the magnetic order in the 
striped-AF phase, vanishing at the transition point. 

The quantities derived from the original order param- 
eters, the cosines of the angle between the neighboring 
spins along the a and b axis, 4> a ^ b \ and the total mag- 
netization |s|, are shown in Fig. [H^b). The behavior of 
cos <fi a and cos <p h confirms the AF order along the b axis 
independent of E z , while in the a direction the angle 
<f) a changes at the phase transition (at E z = 1.27 J) from 
around 27r/3 to n. After the transition the cosines remain 
equal as expected in the isotropic AF phase. The total 
magnetization \s\ is almost constant, increasing monoton- 
ically when E z grows, showing that the essential physics 
of the striped-AF phase lies in the spin angles, although 
its relatively low starting value means that the striped- 
AF phase is affected by strong spin quantum fluctuation 



that weaken AF order. 

In Figs. HPT a) and [TOTb) we present the on-site and 
bond spin-orbital covariances for the same parameter 
range as in Fig. [9j The bond covariance R 1 for a square 
cluster of Fig. (3Ja) is defined as 



i? 7 = ((Si-Si 



lr l T l+7 



>-(S 1 -S 1+7 )(r 1 V7 +7 >, (3 



with 7 = a,b. The striped-AF phase exhibits relatively 
large on-site entanglement in the x spin component, van- 
ishing in the G-AF phase, see Fig. [Hla). In contrast, 
the entanglement in the z spin component persists in the 
G-AF phase, see Fig. |UJb). For the x component the 
dominating covariances are the ones for the b axis and 
for the z component those along the a axis. The fact 
that the z covariances remain finite in the G-AF phase is 
somewhat surprising as one could expect that this phase 
with no frustration and almost fully polarized FO config- 
uration could be trivially factorized into spin and orbital 
wave functions. This expectation based on the previous 
experience^ turns out to be incorrect and we show in Sec. 
IIV CI that high order orbital fluctuation are essential for 
stabilizing the AF order along the c axis in this phase. 



E. Orbital fluctuations 

To estimate the strength of the orbital fluctuations in 
the 3D KK model one can evaluate the total orbital mo- 
ment, defined in a similar way as the total magnetic mo- 
ment of Eq. (J33, i.e., 



1 



z\2 



T x\2 



2V3 



(t a ) 2 + (t b ) 2 +t a t b . 

. (3.9) 

We investigate its value for a representative cut in the 
phase diagram of Fig. [4] taking 77 = 0.13, a realistic value 
for KCuF 3 pi and for — 2J < E z < 1.5J within the cluster 
MF and the single-site MF approximation, see Fig. [TT] 
As expected, for algebraic reasons the orbital moment 
(|3.9[) in the latter approach is trivial - |t| = 0.5 for all 
values of E z . On the contrary, the moment |t| found in 
the cluster MF is reduced from the above maximal clas- 
sical value in all phases except for the ^4-AF phase where 
this reduction is marginal. Quantum phase transitions 
for increasing E z are marked either by discontinuities in 
|r| (first order transitions) or by discontinuities in the 
derivative of |t| (second order transitions). 

The reduction of |r| is most pronounced in the ortho- 
G-AF, where orbital QF couple to spins, and in the PVB- 
AF phase where the continuous orbital phase transition 
takes place but still it does not exceed 20%. We observe 
that the orbital order in the 3D KK model is robust and 
stable against weak QF in all phases. This result follows 
from the rather classical directional nature of e g orbitals 
which leads to the reduction of QF in the orbital spaced 
The orbital order found here in the entire phase diagram 
justifies the pcrturbative expansions in the orbital sector 
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where for simplicity we use a dimensionlcss parameter, 
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FIG. 11. (Color online) Orbital moment |t| (|3~9|) found in the 
cluster MF approximation at the horizontal cut of the phase 
diagram of Fig. [4] for r\ = 0.13 and increasing E, (solid line). 
Quantum fluctuations reduce r from the classical value of 
0.5 found in the single-site MF (dashed line). Quantum phase 
transitions are indicated by vertical dotted lines. 



which are used in Sec. IIVI to derive effective spin mod- 
els. We note that this case is different from t% g orbitals, 
where orbital liquid was found both for the perovskite 
lattice^ and for the frustrated triangular lattice^ for 
the occupancy of one electron per site. 



IV. EFFECTIVE SPIN MODELS 

In this Section we describe effective spin models ex- 
plaining the origin of the exotic magnetic phases. Each 
model is obtained by perturbative expansion around the 
ground state with orbital order stabilized by the domi- 
nant orbital Hamiltonian Ho- The expansion eliminates 
the orbital degrees of freedom yielding an effective spin 
Hamiltonian H s having the ground state with the exotic 
magnetic order. 



A. The ortho-G-AF phase 

We begin with the exotic magnetic order found in the 
ortho-G-AF phase shown in Fig. [5] and derive an effec- 
tive spin model for this phase following the same ideas as 
those employed in the 2D KK model (see Ref. [35h . The 
idea is to use the standard quantum perturbation theory 
with degeneracy for the orbital sector of the KK Hamilto- 
nian. We can divide the Hamiltonian given by Eq. (|2.2[) 
into the unperturbed part Jio and perturbation V in the 
following way: 



H = -Je z ^- 

i 

v = n-n , 



(4.1) 
(4.2) 



J 



(4.3) 



This can serve as a starting point for the perturbative 
treatment for large \E Z \ > J. For negative E z (from now 
on in units of J) the ground state |0) of Ho is the state 
with all z orbitals occupied by the holes, i.e., 



1 



V*: rt |0> = — — |0> , 



(4.4) 



with energy per site £o = ^Je z . Following the quantum 
perturbation theory we can construct the effective spin 
Hamiltonian H s using the expansion in powers of e~ l : 



H s =Ne - 



n^0 



(n\V\ 



f 



-0(e- 



(4.5) 



where all the overlaps are taken between the orbital states 
leaving the spin operators alone, £ n = \E n — Eq\ is the 
excitation energy in the physical units (ex Je z ), and N 
is the number of sites. Knowing the definition of the 
orbital operators Eq. (|2.5[) . we can easily calculate 
the desired orbital overlaps. The first order gives, up to 
the constant term, 



Jo (1) 

j y a b 



JgW (Si • S i+C ) , 



where g^J = (-3ri 



(Si ' Si +7 ) 

i,7— a,b i 

(4.6) 

-4r 2 +r 4 )/2 5 , g£ ] = (r 2 + r 4 ) /2 
are the in-plane and interplane coupling constants. Note 
that <jc is positive in the whole physical range of 77, 
while g^2 is changing sign at rjo ~ 0.1547. If the first 
order term ifi alone were the only spin interaction, then 
•qo would be the point of a quantum phase transition 
between the G-AF and A-AY phases. However, we have 
found that these two phases are separated by a stripe of 
the exotic ortho-G-AF order. 

Since vanishes at 770 and therefore the first order 
in-plane interactions can be arbitrarily weak around 770, 
it is necessary to go to higher order terms of the pertur- 
bative expansion to determine the in-plane interactions 
leading to the exotic magnetic order. In the second order 
the sum runs over all excited orbital states, but from the 
superexchange terms, Eqs. (|2.2p and (|2.3j) . one observes 
that V has non-zero overlap only with states either with 
one or with two orbitals being excited on the considered 
NN bond. This brings us to the second order correction 
of the form: 



H 



(2) 



JfJ 



(2) 



(4.7) 

with gW = 3(r x +2r 2 + 3r 4 ) 2 /2 10 - Here {{ij)) ab and 
(((ij)))ab denote pairs of next NN (NNN) and third NN 
(3NN) sites, respectively, in the a£>-plane (for details of 



11 




FIG. 12. (Color online) Schematic view of the order imposed 
by the Hamiltonian (|4.7p in the ab plane. Interactions be- 
tween NNN spins 1 and 3 (and equivalent) are AF, between 
3NN 1 and 4 are FM, but they vanish for NN spins 1 and 2, 
leaving the angle ip undetermined. 



this deriviation see Appendix |A"]) . The second order ex- 
pansion also yields an ej 1 correction to the NN inter- 
action strength g^V , but this only moves the transition 
point rjQ and does not change the magnetic order. 

One can easily see that the ground state of the Hamil- 
tonian (|4.7|) consists of two antiferromagnets on two in- 
terpenetrating sublattices of the square lattice, stabilized 
by the NNN AF Heisenberg interaction on each sublat- 
tice. Its ground state has long-range AF order renormal- 
ized by large QF. The additional 3NN FM interaction 
makes this AF state more classical and closer to the Neel 
order by suppressing partly QF. In the absence of NN 
coupling at 770, the two AF states on sublattice are uncor- 
rected and the angle ip between the respective AF order 
parameters remains undetermined, as shown in Fig. 1121 
Such ground state is similar to the ortho-G-AF state of 
Fig. [5] with AF interaction along the c axis governed by 
g c , but we still have to find the reason why the magnetic 
moments in two antiferromagnets prefer to be perpendic- 
ular to each other. 

To answer this question we have to go to the third 
order of the perturbative expansion given by 



^i 3) = EE< 



V|n)i(n|V|m)_L(m|V|0) 



(4.8) 

After long but elementary calculations, presented in more 
detail in Appendix [Al and under the assumption that the 
order shown in Fig. [T^] exhibits only weak QF (see Ref. 
|35| for the arguments supporting this statement in the 
2D model based on the spin-wave theory, which apply 
even more here in higher dimension) we obtain a classical 
expression for the energy per site related with the angle 

£i l ] ~ J ^T72 (n + n) (n+ 2r 2 + 3r 4 f cos 2 p ■ (4.9) 
This expression arises from a third-order four-spin inter- 



action. As we can easily see, it is minimized by ip = n/2 
which implies the perpendicularity of the NN spins. This 
finally explains the origin of the ortho-G-AF spin order. 



B. The canted- A- AF phase 

To explain the exotic magnetic order found in the 
canted- A-AF phase, displayed in Fig. [l^b) and described 
in Sec. IIII CI we derive an effective spin model for this 
phase near the orbital degeneracy at E z — in a similar 
way as described above for the the ortho-G-AF phase. 
This time the starting unperturbed orbital Hamiltonian 
stabilizes the AO order, 



U% = Je 



7=0,6 (i,j)\\-f 



E 



(4.10) 



where we gather all the erf cr? terms from the full KK 
Hamiltonian % in Eq. (|2.2|) (plus all constant terms 
which are not relevant and omitted here). This sets the 
coupling constant in Eq. (|4.10|) as 



2 T 



(3n + n) 



(4.11) 



The perturbation is all the rest: V x = H — Hq. The 
ground state of r hi% is a classical configuration with the 
alternating eigenstates erf. This order is consistent with 
the cluster MF results showing that the canted- A-AF is 
an AO phase, just like the neighboring A-AF and FM 
phases. 

The effective spin Hamiltonian H s can be constructed 
using the formal expansion of Eq. (|4.5[) . In the zeroth 
order we get the ground state energy of Hq, and in the 
first order the spin Hamiltonian: 

Hi 1] =-J 9 ^ E (S i -S J+7 ) + J.gW^(S J .S J+c ), 

i,y= a, b i 

(4.12) 

with S £= 7(-8r 2 + 7r 1 -r 4 )/2 5 , = (2r 2 -n+r 4 )/2 3 . 

The NN interaction along the c axis g£\rj) is changing 
sign at j)l ~ 0.236. This agrees well with the AF-FM 
crossover in the c axis taking place in the canted- A-AF 
phase that separates the A-AF and FM phases for 77 close 
to ?7q. The NN interaction in the ab planes g^ is FM 
in this region. Since the first order interaction along the 
c axis can be made arbitrarily weak near t]q, we have 
to go to higher orders to describe interactions between 
different ab planes. 

In the second order V can produce two types of ex- 
cited states: (i) with a single orbital rotated by 7r/2 in 
the ab plane, or (ii) with a rotated pair of neighboring or- 
bitals. This leads to two type of terms in the interplanar 
Hamiltonian, 



(2) 



Jg, 



(2) 



8e* 



8e* 



(4.13) 
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with g^ = {§ (r? - rj) - (r 2 + r 4 ) 2 }/2 6 + ^(r, +r 4 ), 
and gcc = (f2 + r 4 ) 2 /2 7 (for more details see Appendix 
|B1) . The NN Hciscnberg interaction <?c renormalizes 
again the first order term modifying slightly the 

value of t]q where the NN coupling changes sign and 
makes it dependent on E z . The ^-dependence of the 
crossover from the A-AF to FM phase is consistent with 
the phase diagram in Fig. 2] The NNN interaction gic is 
FM in agreement with the FM NNN correlations in the 
canted magnetic order shown in Fig. |6|b) and supported 
by the numerical result in Fig. [5Ja) . The nontrivial cant- 
ing angle 9, interpolating between 9 = n and 9 = when 
passing from the A-AF to FM phase, remains undeter- 
mined and we must proceed to the third order of the 
perturbative expansion, see Appendix [Bj 

The perturbative expansion up to third order leads to 
the classical expression for the ground state energy: 



E (9) _ cos 9 

j ~ ~T~ 

and we write 



,(!)- 



8c 



-0 



cos 2 9 1 
32 T 2 / 



; (3) 



E (9) = JA(rj) cos 9 + JB(rj) cos 2 9 . 



(4.14) 



(4.15) 



Since B(ri) > 0, the energy has a nontrivial minimum at 
9 given by the equation: 

( -1 for A(rj) < 

cos9 Q =l-2^ for \A( V )\ < fflfa) (4.16) 
{ 1 for A(rj) > \B{n). 

This reproduces well the behavior of cos 9 obtained via 
the cluster MF method, see Fig. [TJa). 



C. The G-AF versus G-AF spin order 

Here we derive an effective spin Hamiltonian around 
the FO ordered state with x orbitals occupied by holes 
to explain the energy difference between the G-AF and 
C-AF phase, favoring the G-AF order in the cluster MF 
approach, see Fig. |U As in Sec. IIV Al we can divide 
the Hamiltonian given by Eq. (|2.2[) into the unperturbed 
part Ho (|4.1[) and the perturbation V consists of intersite 
terms in the ab planes (V a b) an d along the c axis (V c ), 



V = J(V ab + V c ) 



(4.17) 



They are given as follows: 
V ab = -\ £ H7 i+J , V c = ~J2 H ti+o- ( 4 - 18 ) 



z,7— a,b 



For positive E z the ground state |0) of Ho is the state 
with all x orbitals occupied, i.e., 



V*: 7f|0) = -|0), 



(4.19) 




FIG. 13. (Color online) Schematic view of the G-AF spin 
order on two ab planes of 3D cubic lattice stabilized by the 
third-order FM NNN Heisenberg interactions between sites 
i + j and i + c shown by dashed (green) lines (see Eq. (|4,22[) ~). 
The NN FM interactions between i and i+c, shown by vertical 
(red) line, are frustrated in the G-AF phase. 



with energy E = —\Je z per site. Using Eq. (|4.5[) we 
construct the effective spin Hamiltonian H s as a power 
expansion in V. Note that the ground state |0) is an 
cigenstatc of V c to zero eigenvalue; for this reason V c 
gives no contribution to H s up to second order in V and 
the interplane interactions appear as e~ 2 terms. 

The effective spin Hamiltonian H s up to second order 
in V can be calculated in the same way as in Sec. IIV Al 
and contains the same spin interactions, 



H. = j{ g m+o( c ; 1 )} VJ(s,-s, 



<ij>„6 ««». 

^ J2 (s.-s.o + o^ 2 ), 



(4.20) 



fJ 



(i) 



3(4r 2 - n + 3r 4 )/2 5 , = 3(2r 2 + r 4 - r 1 ) 2 /2 w 1 
and where all bonds are in the ab planes. Again, one can 
think of an ortho-G-AF phase similar to the one obtained 
in the limit of E z — > — oo for rj close to rjo where the NN 
interaction g^ vanishes (in this case tjo — 0.2915) but 
this is not our aim here — we search for an effective 
description of the interplane interactions. 

Looking at the third order correction to H s given in 
Eq. (|4.8|) and at the form of V c we can see that the 

(3) 

interplane part of the third order correction HsJ has the 
following structure: 



H (3) 

s,c 



\Vab\n) ^ {n\V c \n) — (n|V a6 |0) 



(4.21) 

This can be greatly simplified if we take into account 
two facts: (i) V a b can excite only single orbitals or pair of 
neighboring orbitals so the sum over n turns into the sum 
over excited sites i, and (ii) V c is non-zero only around the 

excited orbitals. The final form of i?i 3 c can be obtained 
using some identities for spin products and assuming that 
the spins form classical AF state in the ab planes (see 
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Appendix [C]). Of course, this is not exact because we 
neglect QF but these should be small in the 3D long- 
range ordered magnetic configuration being either G-AF 
or C-AF. The final formula for interplane interactions is: 



Hi 3 ) = -J 



7 



, E {^ 3) (* 

>}• 



(S 



i-\-c 



(4.22) 



with g { c ] = n (n - r 4 ) and = i (n + r 4 ) (3n - r 4 ) 
containing NN and NNN Heisenbcrg bonds both favoring 
fcrromagnctism (S r takes only non- negative values). 

This could lead to frustration as the NN bonds favor 
the C-AF state while the NNN bonds prefer the G-AF 
one, but one can easily check that the number of NNN 
bonds per site is four times larger than that of NN bonds 
(see Fig. [i~3|) . If we take into account this factor of 4 
comparing the coupling constants gi 3 ^ and g^J of the 
NN and NNN interactions we obtain: 



9, 



(3) 



1 + 77 (1 - 4??) 
(3?7 2 + 2 ?? - l) 2 



> 



(4.23) 



in the physical range of 77 £ [0, 1/3). This finally explains 
why the G-AF phase is favored over the G-AF phase, as 
long as we go beyond the single-site MF approximation 
which cannot capture the subtle third-order orbital fluc- 
tuations (compare Figs. [2] and [4]). 



D. The striped- AF phase 



In Sec. IIIIDI we summarized cluster MF results for 
the striped- AF phase and the continuous phase transition 
between the striped- AF and G-AF phases. These results 
show in particular that, unlike all other exotic phases, 
the striped phase has relatively large QF in both the 
orbital and the spin sector. Nevertheless, we attempt 
explaining its origin by an effective spin model keeping 
in mind qualitative character of our analysis. 

The perturbativc expansion follows the same lines as 
in Sec. IIV CI We assume that E z > and the unper- 
turbed ground state |0) is fully polarized by E z , see Eq. 
(|4.19p . This is not an unreasonable starting approxima- 
tion, because in Fig. Eta) we find t c = (rf) rj 0.38 in the 
striped phase, i.e., only about 12% of orbitals are flipped 
with respect to the fully polarized state |0). We believe 
that it is still justified to make a perturbativc expansion 
in the orbital sector around the FOx state. We proceed 
with the expansion in the same way as in Sec. but 
here we assume AF order along the c axis and focus on 
the ab planes only. 

Up to second order the effective Hamiltonian is given 
by Eq. (|4.20[) which can be rewritten in a more compact 
form, 



H s — Ji 



E 



I (Si-s,o-i e &- s A 

[<«» (««»> J 



,(4.24) 



with both Ji , J2 > in the physically interesting range 
of 77 and E z . The first two terms are the same as in the 
Ji — J2 model whose ground state is AF when J\ ^s> J2 
and collinear when J 4 <g; J2 ■ The extra 3NN FM coupling 
is consistent with the NNN AF coupling, hence the phase 
diagram of the present model should be qualitatively the 
same. 

In a classical approximation, the phase diagram of the 
Ji — J2 model (|4.24p motivates an Ansatz, where 



S,;-S 



i+b 



1 



Si-S i+a = ^cos0 a = -icosa, (4.25) 

with a = it — <fi a , sec Fig. EJa). Here the angle a is a 
variational parameter. The Ansatz breaks the symmetry 
between the a and b axes. It can describe the AF phase 
when a = 0, the collinear phase when a = 7T, and the 
intermediate striped- AF phase when a € (0,7r). Up to 
an additive constant, the energy per site is 



•(a) = (-,/!+ 2J 2 ) 



cos a . 



(4.26) 



It suggests a discontinuous (first order) transition be- 
tween the AF and collinear phases at J\ = 2J 2 . Thus, 
up to second order in the pcrturbative expansion, the 
striped phase is unstable in the classical version of the 
effective spin Hamiltonian H s . However, near the criti- 
cal point, where the energy e(a) does not depend on a, 
the ground state can be very susceptible to any pertur- 
bation from higher order terms in H s . 

Indeed, the third order term Hs contains the relevant 
perturbation (|4.8|) of the form 



H s%b — ^3 E ( S » ' S i) E S 1' S 1" ( S i+7' ' S J+7") ' 

(4.27) 

where 7', 7" stand for possible direction {a, 6} in the 
square lattice, and s± a r = 1, s±b = —1, respectively. 
Including these interactions, the energy per site becomes 

e(a) = (-Ji + 2 J 2 + 4 J 3 ) cos a + 5J 3 cos 2 a. (4.28) 

Now the system undergoes a continuous symmetry- 
breaking phase transition from the AF phase (\a\ = 0) 
to the striped phase (|a| > 0) when Ji becomes less than 
2J2 + 14 J3. This is where the quadratic term 62 in the 
Landau expansion of energy, s(a) pa £0 + £2« 2 + 84a 4 , 
becomes negative favoring a finite value of the order pa- 
rameter a. The continuous character of the transition is 
consistent with the cluster MF results in Fig. |H1 

In conclusion, the effective classical spin model pro- 
vides a qualitative explanation for the origin of the 
striped-AF phase. We do not present any quantitative 
predictions here — they would be rather poor due to 
large spin QF in the striped phase. 
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V. SUMMARY AND CONCLUSIONS 

We have presented a rather complete analysis of the 
phase diagram of the 3D Kugel-Khomskii model in the 
framework of the cluster MF theory and effective pertur- 
bative models. It is found that spin disordered plaquette 
valence-bond phase is stable near the orbital degener- 
acy in the broad range of Hund's exchange, resolving the 
existing controversy and in agreement with the earlier 
studicsiii2S Furthermore, we managed to obtain a very 
transparent picture in the left part of the phase diagram 
of Fig. 2] (for negative crystal- field splitting) , with a se- 
quence of phase transitions occurring for increasing r\ and 
involving two intermediate phases with exotic magnetic 
orders. This sequence is caused by the (first order) NN 
Heisenberg interactions changing their sign from AF to 
FM spin interaction: (i) first in the ab plane, where the 
ortho-G-AF phase occurs, and (ii) next along the c axis, 
which leads to the canted- A-AF configuration. In both 
cases we found perturbativc expansions around certain 
orbital configurations explaining these puzzling magnetic 
orders by further neighbor spin interactions. In case of 
the ortho-G-AF phase the effective Hamiltonian is the 
same as for the 2D KK model with additional AF interac- 
tions along the c axis. In the other case the intermediate 
configuration turned out to be essentially classical, with 
FM planes damping spin quantum fluctuations, and we 
have used this fact to construct the effective spin model 
around the AO configuration to explain the canted- A- AF 
ordering. 

To supplement these analytical considerations we pre- 
sented the plots of order parameters, correlations and 
spin-orbital covariances for the two cuts in the phase 
diagram, passing through canted- A- AF and striped- AF 
phases. They show that both phases involve spin-orbital 
entanglement^ and the plots for the canted- A- AF phase 
are supported by the effective spin Hamiltonian derived 
for this phase. The plots for the ortho-G-AF phase were 
already given in the 2D case,— and they are rather sim- 
ilar for the present 3D cubic lattice (not shown). 

In contrast, some of our results in the right part of 
the phase diagram (for positive crystal-field splitting) 
arc more qualitative. As the last analytical result we 
showed the derivation of the effective spin Hamiltonian 
for the FOa; configuration which explains why the G-AF 
phase always wins over the G-AF order for the 3D and 
bilayer systems although these two phases are degenerate 
in the single-site MF approach. The answer was found 
in the third order of the perturbation expansion and it 
was showed that the AF order along the c axis is induced 
by the two effects: (i) AF order in the ab planes, and 
(ii) NNN interaction between the planes being FM. Sur- 
prisingly, this demonstrates that the right G-AF phase 
is highly spin-orbital entangled, as shown in Fig. [TU] be- 
cause we need third order orbital fluctuations to stabilize 
the spin order, and agrees with the plots of spin-orbital 
covariances presented for the bilayer KK model^ This 
may be the reason for a stronger divergence of the quan- 



tum corrections to the order parameter found in the spin- 
wave theory^ 

Our study has shown that the striped- AF phase still 
requires a more sophisticated approach than the one that 
was implemented here. In particular, it is not clear why 
this phase appears only in the 3D case; we have verified 
that it converges to a stable solution in the 2D case, but 
with energy being always higher than that of either the 
G-AF or FM phase. This remains one of the open ques- 
tions in the phase diagram of the 3D Kugel-Khomskii 
model and further studies beyond the cluster MF, such as 
the entanglement renormalization ansatz^ are needed. 

Summarizing, we have established that orbital exci- 
tations not only couple to spin fluctuations^ but may 
even change the spin order in spin-orbital systems. On 
the example of the 3D Kugel-Khomskii model we have 
shown that three magnetic phases with exotic spin or- 
der arise near the crossover from AF to FM spin inter- 
actions at increasing Hund's exchange and are triggered 
by entangled spin-orbital quantum fluctuations. The de- 
rived effective spin models, that include second and third 
neighbor spin interactions and go beyond the Heisenberg 
paradigm, provide good microscopic insights into these 
phases. At the same time, these effective interactions 
do not confirm the earlier suggestion^ that the resonat- 
ing valence-bond phase accompanied by the FO order of 
3z 2 — r 2 orbitals along the c axis could be stable instead. 

The derived phase diagram could in principle de- 
scribe also the A-AF phase observed in KCUF3 at low 
temperature jil taking realistic parameters, but the sig- 
natures of the ID Heisenberg physics are missing. On 
the one hand, it supports the recent view that the Kugel- 
Khomskii model is incomplete and should be extended by 
the Goodenough processes^ and by other (lattice) degrees 
of freedom^ to explain fully the observed physical prop- 
erties of KCUF3. On the other hand, the present study 
provides an experimental challenge whether the ortho- 
G-AF phase could be discovered in KCUF3, for instance 
using high pressure studies. 
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Appendix A: Spin model in the ortho-G-AF phase 

The form of second order contribution to the effective 
in-plane spin Hamiltonian can be derived from Eq. (|4.5|) 
by calculating the matrix elements (n| V |0) in the orbital 
sector. As a result we get 

^=^ (2) e{e^(s«-s j+7 )| 
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FIG. 14. (Color online) Schematic views of second order cor- 
rections in the effective spin Hamiltonian: (a) effective NNN 
interactions, and (b) effective 3NN interaction. Red frames 
stand for Heisenberg bond with ± sign depending on the 
bond's direction and filled (magenta) dots are the single-site 
orbital excitations in the ground state. 



w) 



ri + r 4 ) + 



3ri - r 4 



(Al) 



where s 7 is a sign factor depending on the bond's direc- 
tion 7 and originating from the definition of operators 



a(b) 

T- 1 9. 

I , i.e., 



1 if 7 = ±a 
-1 if 7 = ±fe 



(A2) 



and is defined in Section IIV Al The squared quan- 

(2) 

titles in Hs produce spin products of the two forms, 
shown in Fig. fT4Ta) and ITW b) . which can be simplified 
using elementary spin identities: 



(Si-f^ ■ S{) (Si ■ Sj-py') — (Sj- 



-7 ■ Sj_|- 7 ') 



'i+7 



(Si X Sj+y) ,(A3) 



(Si • Si+ 7 ) 



-I ( Si • S i+T ) + A . (A4) 

The latter identity simplifies the second line of W ' ab and 
produces E~ x correction to the Heisenberg interactions 

in Hs , while the former one leads to the interactions 
between further neighbors in Eq. (|4.7j) . The imaginary 
term in Eq. (|A3|) is antihermitian and must cancel out 



(2) . (2) 

with other terms in H^ ab . Thus, in ab , we arc left with 

the pure Heisenberg term (Si+ 7 • Si+y) connecting sites 

% + 7 and i + 7', being either NNN or 3NN, presented 

in Fig. IT4T a) and [F4T b). with the sign given by s 7 s 7 '. 

Due to the double counting of the interactions the 3NN 

couplings in H^ ab of Eq. (|4.7| is twice weaker than the 

NNN ones. 

The third order correction necessary to determine the 
in-plane NN interaction at 770 is given by Eq. (|4.8[) ; 
here we limit the sums over the excited orbital states 
{|n) , \m}} only to the ones with orbital flips lying in 
the same ab plane. This produces many contributions 




FIG. 15. (Color online) Schematic view of third order cor- 
rections to H 3 (frames stand for Heisenberg bonds with ± 
sign depending on their direction, filled circles indicate or- 
bital flips): (a) term bringing new physics, chain product of 
three bonds, (b) term doubling the results from second order 
expansion, interaction between spins 2 and 3 is squared (bond 
marked by a thick frame) . 



to the spin Hamiltonian but we are interested only in 
terms with new operator structure with respect to lower 
orders because the others will be just the E~ 2 corrections 
to the already existing interactions. The terms bringing 
potentially new physics are the ones with three different 
Heisenberg bonds multiplied one after another. Such con- 
tribution is depicted in Fig. fTSTa) for sites i = 1,2,3,4 
(in contrast, Fig. [TSTb) shows the term bringing no new 
contribution to H s ) and can be transformed following the 
identity, 

(S r S 2 )(S 2 -S 3 )(S3-S 4 ) = 

ls 1 -S4 + ^(S 1 -S 4 )(S 2 -S 3 ) 
16 4 

-i(S 1 -S 3 )(S 2 -S 4 ) + gS 1 -(S3 x S 4 ) 

+ ^S 1 .(S 4 xS 2 ) + ^S 1 .(S 4 xS 2 ), (A5) 
o o 

where the cross-product terms are antihermitian and 
must cancel out with other terms of the same structure 
in H^lb- To analyze the remaining terms in Eq. (|A5[) 
it is helpful to employ the almost classical nature of the 
AF spin order on two sublattices: (i) the first term is an 
AF interaction between the sublattices which is not com- 
patible with the antiferromagnetism on sublattices that 
is one order of E z stronger, (ii) the second term favors 
perpendicularity of the two AF sublattices, as long as its 
sign is positive, which is compatible with the order on 
sublattices, and (iii) the third term brings no new infor- 
mation about the spin order. 

Putting all together, we can argue that the third order 
perturbative contributions in Eq. (|A5|) can favor per- 
pendicularity of the order parameters on two sublatticc 
antiferromagnets. Now we have to extract all such contri- 
butions from Eq. (|4.8|) and check whether the total sign is 
indeed positive because otherwise our arguments would 
be incomplete. After lengthy but elementary calculation 
(see Ref. for more details) we obtain the third or- 
der Hamiltonian, with interactions not included in lower 
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orders of the form: 

H i% = J TO ^ + r *) ^ + 2r 2 + 3r 4 ) 2 



2 17 e 

E 

<*J")II7 



(Si • Sj) S ^ Sj'Sj" (Sf_py/ *Sj-py") . 



7 7*7 



(A6) 



These are the special contributions of the type shown in 
Fig. rT5T a) taken into account here. This, in the clas- 

sical limit, gives the energy £j_ in Eq. (|4.9[) favoring 
perpendicular orientation of the NN spins. 



Appendix B: Spin model in the canted- A-AF phase 

To calculate the second order correction to the spin 
Hamiltonian H s in the canted- phase we have to, as 
before, calculate the matrix element (overlap) (n\ V x |0) 
in the orbital sector. Since the state |n) is an excited 
state of Hq, the non-zero overlap can be obtained only 
for the part V x of the interaction V that contains at least 
one operator erf. Assuming FM order in the ab planes it 
can be written as, 



E {^K-s 7 V3(crK + ^K)} 

m\\ b 

<«>u 

(r 2 + r 4 )H^> (at + of) } - \ E z £ a\ . (Bl) 



employing the spin projectors in Eqs. (|2.6j) . Now we can 
derive the second order Hamiltonian, 



H 



(2) 




3e x 



E^H 



3ri - r 4 



(B2) 



which can be easily transformed into Eq. (|4.13[) using 
the identity (|A3)| . 

In the third order we are interested only in terms with 
three Heiscnbcrg bonds multiplied by each other along 
the c direction, according to Eq. (|A5|) . as all other pos- 
sible products will only reproduce the results from lower 

(3) 

orders. Such terms lead to the correction HsJ of the 
form: 

t„(3) 

Hi 3 ) = V{(S J+C • Si) (Sj • S,_ c ) (S 4 _ c • Si_ 2c ) + H.c.} , 



r2 1 



(B3) 



(3) _ 5 /_ , _ n2 
? c ~ 6 

ther simplihed using the identity (|A4[) . 



with ,^ J; = | (r 2 + r 4 ) z (r 4 + n) /2 14 . This can be fur- 



(3) 

x I 

-4(S i+c -Si_ c )(Si-Si_ 2c )}. (B4) 

Now, according to Fig. |6jb) , we use the classical expres- 
sions for the spin scalar product using the canting angle 9 
for the odd neighbors and FM order for the even neigh- 
bors, imposed by the second order Hamiltonian of Eq. 
6J31), i.e., 

Si • S i+ (2„_i) c = — COS (9, Si ■ Si-|_2nc = — . (B5) 



After inserting this into Eq. (|B4|) only the first line 
depends on 8 and gives a contribution to the classical 
ground state energy Eo(9) of the canted- ^4-AF phase see 
Eq. f£I5|) . 



Appendix C: Spin model in the G-AF phase 

Here we consider the G-AF phase with FOx order for 
E z > 0. The third order contribution of the form given 
by Eq. (|4.21[) can be expressed as, 



Hi s 3 ^ - j9 ~^f E S 7V (Si-S i+7 ) (S r S i+c ) (Si-S 



2s 



^EU^+^)^- s 



2 13 e? 



3r a - r 4 



x (Si-Si+c) ^ (n + r 4 ) (Si-Si +7 ) + 3ri . Ti }, (CI) 



where S r = ?'i — r 2 , sums over 7 and 7' run over the 
in-plane directions {±a,±6}. The first term in Eq. (|C1[) 
comes from single-orbital excitations, while the second 
one (second and third line) — from two-orbital excita- 
tions. The interplane interactions are typically sand- 
wiched between two in-plane bonds; these can be treated 
with the following spin identity: 

(Sj-Si4- 7 ) (Si-Sj-|- C ) (Si-Sj4- 7 ') 

= 4 (Sj+c - Sj+ 7 ) (Sj'Si+y) — — (Sj+ 7 < • Si+ 7 ) (Si-Si-)_ c ) 

1 i 
+ — (Si-Sj+ 7 ) (Si +c -S; +7 ') + — Si +7 -(Si +c x Sj+y) . 
4 o 

(C2) 

Under the assumption of the classical AF order in ab 
planes, i.e., (Si • Si+ 7 ) = — j for 7 = ±a, ±6, the above 
identity gives the following results for the tri-quadratic 
spin products in H^ 3 }, 

X^ 7 , 7 " S 7 S 7" (Si • Si-f 7 ) (Si • Sj-|_ c ) (Sj • Si_)_ 7 ") 

= j E ( s i+c-Si +7 ) + - (Si-Si +C ) , (C3) 
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and 

(Sj-Sj+ 7 ) (Sj-Sj+c) (Sj-Si_)_ 7 ) = -— (Sj-Si+ C ) . (C4) 

The biquadratic terms can be simplified using identity 
(1X3)) . Note that the classical AF order, (S, • S i+7 ) = —\ 



for 7 = ±a, ±6, is inserted after the spin products are 
sorted out with spin identities (|A3[) and (|C2[) but not 
before. This is an important distinction: we extract the 
result of higher order in-plane spin fluctuations before 
we freeze them out to get a simplified picture. After 

(3) 

gathering all the interactions together we obtain HsJ in 
Eq. KTIh . 
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